######################### EXCESS DEMAND FUNCTION

if params.mobility_spec=="baseline" 

    #Excess demand function
    function excessdemand(p_ic, s_ij_c)

        #County commodity prices
        p_i_c_pretax = reshape(p_ic,params.I,params.C)
        p_i_c = p_i_c_pretax - params_policy.tax_i
        AN_i = params.AN_i + params_policy.subsidy_i 

        #Intermediate input prices
        wM_i = params.wM_i
        wM_i_c = repeat(wM_i,1,params.C)

        ################################ SUPPLY
        
        #### Solve for labor market equilibrium
        function excessdemand_H(wH_i)

            wH_i_c = repeat(wH_i,1,params.C)

            if (params_policy.policy=="lf") && (params.back_out_zeta=="y")
                zeta_i_c = params.zetat_i_c.*(wH_i_c.^(params.gamma_H_i./params.gamma_L_i)).*(wM_i_c.^(params.gamma_M_i./params.gamma_L_i))
            else
                zeta_i_c = params.zeta_i_c
            end

            #Labor demand
            A_i_c = params.z_i_c.*zeta_i_c
            r_i_c_p = p_i_c.^(ones(params.I,params.C)./params.gamma_L_i)
            r_i_c_H = wH_i_c.^(-params.gamma_H_i./params.gamma_L_i)
            r_i_c_M = wM_i_c.^(-params.gamma_M_i./params.gamma_L_i)
            r_i_c = A_i_c.*r_i_c_p.*r_i_c_H.*r_i_c_M
            rtilde_i_c = r_i_c./(p_i_c.*params.gamma_L_i)

            sumexp_C = sum(r_i_c.^params.theta_lambda, dims=2)
            pi_cC = (r_i_c.^params.theta_lambda)./sumexp_C
            pi_N = (AN_i.^params.theta)./(AN_i.^params.theta+sumexp_C.^params.lambda)
            pi_C = ones(params.I,1) - pi_N
            pi_c = pi_cC.*pi_C
            L_i_c = pi_c.*params.L_i
            QS_i_c = rtilde_i_c.*(pi_cC.^(-params.theta_lambda.^(-1))).*(pi_C.^(-params.theta.^(-1))).*L_i_c
            H_i_c = params.gamma_H_i.*p_i_c.*QS_i_c./wH_i_c
            M_i_c = params.gamma_M_i.*p_i_c.*QS_i_c./wM_i_c
            H_i_demand = sum(H_i_c,dims=2)

            #Labor supply
            share_i = params.Hbar_i./sum(params.Hbar_i)
            shareHcond_i = (wH_i.^(params.psi))./(wH_i.^(params.psi) + params.wH_i_outside.^(params.psi))
            shareH_i = shareHcond_i.*share_i
            H_i_supply = params.Hbar_i.*( shareHcond_i.^((params.psi-1)/params.psi))

            #Excess demand
            xd_H_i = (H_i_demand-H_i_supply)./(minimum(hcat(H_i_demand,H_i_supply), dims=2))

            return (xd_H_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c)

        end

        #Inner labor market loop begins here
        wH_i = ones(params.I)
        x_a = wH_i.*params_solver.x_a_H
        x_b = wH_i.*params_solver.x_b_H
        iter_H = 1
        iter_H_max = params_solver.iter_H_max
        xd_tol_H = params_solver.xd_tol_H
        xd_H = xd_tol_H + 1
        xd_H_state = "declining"

        while (xd_H > xd_tol_H) && (iter_H < iter_H_max)  && (xd_H_state == "declining")

            f_a, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_a)
            f_b, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_b)

            x_c = (x_a + x_b)/2
            f_c, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_c)

            if (xd_H != maximum(abs.(f_c)))
                xd_H_state = "declining"
            else
                xd_H_state = "stalled"
            end

            xd_H = maximum(abs.(f_c))

            fc_samesign_a = abs.(sign.(f_c)+sign.(f_a))/2
            fc_samesign_b = abs.(sign.(f_c)+sign.(f_b))/2
            
            x_a = x_c.*fc_samesign_a + x_a.*(ones(params.I,1)-fc_samesign_a)
            x_b = x_c.*fc_samesign_b + x_b.*(ones(params.I,1)-fc_samesign_b)

            iter_H = iter_H + 1

            wH_i = x_c

        end

        if (xd_H <= xd_tol_H) && (iter_H <= iter_H_max)  && (xd_H_state == "declining")
            labor_market_cleared = 1
        else
            labor_market_cleared = 0
        end

        if (iter_H == iter_H_max)
            xd_H_state = "maxed-out"
        else
        end


        #### Evaluate at equilibrium wages
        xd_H_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c  = excessdemand_H(wH_i)
        xd_H = maximum(abs.(xd_H_i))


        #### Supply elasticities
        dQdP_aux = (params.theta_lambda-ones(params.I,1)).*(ones(params.I,params.C)-pi_cC) + (params.theta-ones(params.I,1)).*pi_cC.*(ones(params.I,1)-pi_C)
        dQdP = (dQdP_aux./params.gamma_L_i) + (ones(params.I,params.C)./params.gamma_L_i)-ones(params.I,params.C)
        dPdQ= dQdP.^(-1)


        ################ Intermediaries

        #Marginal product
        mp_ij_c = reshape(params.mp,params.I*params.C)*ones(1,params.J)

        #Markdown
        if params_policy.conduct=="pc"
            mu_ij_c = ones(params.I*params.C,params.J)
        else
            s_ij_c = s_ij_c./s_ij_c
            N_ij_c = reshape(params.N_firms,params.I*params.C)*ones(1,params.J)
            mu_ij_c = ( ones(params.I*params.C,params.J) + reshape(dPdQ,params.I*params.C,1).*(s_ij_c./N_ij_c) ).^(-1)
        end


        ################ Demand

        #County-destination-commodity prices
        p_ij_c_pretax = (reshape(p_i_c_pretax,params.I*params.C,1).*repeat(params.tau_ij,params.C,1))./(mu_ij_c.*mp_ij_c)
        p_ij_c = p_ij_c_pretax + params_policy.tax_ij_c

        #Nation-destination-commodity price index
        aux_n = sum(repeat(params.a_ij_c.*(p_ij_c.^(ones(1,params.J)-params.eta_l')),1,params.N*params.C).*kronecker( Diagonal(ones(2)),kronecker(params.indicator_in,ones(1,params.J))),dims=1).^((ones(1,params.J*params.N*params.C)-kronecker(ones(1,params.C+params.N),params.eta_l')).^(-1))
        P_nj_c = reshape(aux_n,params.J,params.N*params.C)'

        #Destination-commodity price index
        aux_j1 = params.a_nj_c.*(P_nj_c.^(ones(1,params.J)-params.eta_m'))
        aux_j = hcat(aux_j1[1:params.N,:],aux_j1[params.N+1:2*params.N,:]) #horizontally concatenate C times
        P_j_c = (reshape(sum(aux_j, dims=1),params.J,params.C).^((ones(params.J,1)-params.eta_m).^(-1)))'
        
        #Destination price index
        P_j = sum(params.a_j_c.*(P_j_c.^(ones(1,params.J)-params.eta_u')),dims=1).^((ones(1,params.J)-params.eta_u').^(-1))

        #County-level demand
        aux_ij_c = params.a_ij_c.*(p_ij_c.^(-params.eta_l'))
        aux_nj_c = (params.a_nj_c.*(P_nj_c.^(params.eta_l'-params.eta_m')).*kronecker(params.a_j_c.*(P_j_c.^(params.eta_m'-params.eta_u')),ones(params.N,1))).*(P_j.^(params.eta_u'-ones(1,params.J))).*params.X_j'

        C_ij_c = aux_ij_c.*(kronecker(Diagonal(ones(params.C)),params.indicator_in)*aux_nj_c)
        Q_ij_c = C_ij_c.*repeat(params.tau_ij, params.C)
        QD_i_c = reshape(sum(Q_ij_c, dims=2),params.I,params.C)
        
        #Source shares
        s_ij_c = Q_ij_c./sum(Q_ij_c, dims=2)

        #Excess demand
        ED_ic = reshape(QD_i_c - QS_i_c,params.I*params.C)
        QDQS_den_max = maximum( hcat(reshape(QD_i_c,params.I*params.C),reshape(QS_i_c,params.I*params.C)), dims=2)
        QDQS_den_min = minimum( hcat(reshape(QD_i_c,params.I*params.C),reshape(QS_i_c,params.I*params.C)), dims=2)

        xd_ic_max = ED_ic./QDQS_den_max
        xd_ic_min = ED_ic./QDQS_den_min
        
        return (xd_ic_max, xd_ic_min, pi_cC, pi_C, pi_c, pi_N, p_i_c_pretax, p_i_c, p_ij_c_pretax, p_ij_c, QD_i_c, QS_i_c, C_ij_c, Q_ij_c, s_ij_c, P_j_c, sumexp_C, dQdP, dPdQ, mu_ij_c, labor_market_cleared, iter_H, wH_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, xd_H, xd_H_i, xd_H_state)

    end

else
end

if params.mobility_spec=="mobility"

    #Excess demand function
    function excessdemand(p_ic, s_ij_c)

        #County commodity prices
        p_i_c_pretax = reshape(p_ic,params.I,params.C)
        p_i_c = p_i_c_pretax - params_policy.tax_i
        AN_i = params.AN_i + params_policy.subsidy_i 

        #Intermediate input prices
        wM_i = params.wM_i
        wM_i_c = repeat(wM_i,1,params.C)

        ################################ SUPPLY
        
        #### Solve for labor market equilibrium
        function excessdemand_H(wH_i)

            wH_i_c = repeat(wH_i,1,params.C)

            if (params_policy.policy=="lf") && (params.back_out_zeta=="y")
                zeta_i_c = params.zetat_i_c.*(wH_i_c.^(params.gamma_H_i./params.gamma_L_i)).*(wM_i_c.^(params.gamma_M_i./params.gamma_L_i))
            else
                zeta_i_c = params.zeta_i_c
            end

            #Labor demand
            A_i_c = params.z_i_c.*zeta_i_c
            r_i_c_p = p_i_c.^(ones(params.I,params.C)./params.gamma_L_i)
            r_i_c_H = wH_i_c.^(-params.gamma_H_i./params.gamma_L_i)
            r_i_c_M = wM_i_c.^(-params.gamma_M_i./params.gamma_L_i)
            r_i_c = A_i_c.*r_i_c_p.*r_i_c_H.*r_i_c_M
            rtilde_i_c = r_i_c./(p_i_c.*params.gamma_L_i)

            sumexp_C = sum(r_i_c.^params.theta_lambda, dims=2)
            pi_cC = (r_i_c.^params.theta_lambda)./sumexp_C
            pi_N = (AN_i.^params.theta)./(AN_i.^params.theta+sumexp_C.^params.lambda)
            pi_C = ones(params.I,1) - pi_N
            pi_c = pi_cC.*pi_C
            L_i_c = pi_c.*params.L_i
            QS_i_c = rtilde_i_c.*(pi_cC.^(-params.theta_lambda.^(-1))).*(pi_C.^(-params.theta.^(-1))).*L_i_c
            H_i_c = params.gamma_H_i.*p_i_c.*QS_i_c./wH_i_c
            M_i_c = params.gamma_M_i.*p_i_c.*QS_i_c./wM_i_c
            H_i_demand = sum(H_i_c,dims=2)

            #Labor supply
            V_i = params.nu_i.*((wH_i.^(params.psi) + params.wH_i_outside.^(params.psi)).^(1/params.psi))
            share_i =  (V_i.^params.varphi)./(sum(V_i.^params.varphi))
            shareHcond_i =  (wH_i.^(params.psi))./(wH_i.^(params.psi) + params.wH_i_outside.^(params.psi))
            shareH_i = shareHcond_i.*share_i
            H_i_supply = sum(params.Hbar_i).*(shareH_i.^((params.psi-1)/params.psi))

            #Excess demand
            xd_H_i = (H_i_demand-H_i_supply)./(minimum(hcat(H_i_demand,H_i_supply), dims=2))

            return (xd_H_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c)

        end

        #Inner labor market loop begins here
        wH_i = ones(params.I)
        x_a = wH_i.*params_solver.x_a_H
        x_b = wH_i.*params_solver.x_b_H
        iter_H = 1
        iter_H_max = params_solver.iter_H_max
        xd_tol_H = params_solver.xd_tol_H
        xd_H = xd_tol_H + 1
        xd_H_state = "declining"

        while (xd_H > xd_tol_H) && (iter_H < iter_H_max)  && (xd_H_state == "declining")

            f_a, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_a)
            f_b, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_b)

            x_c = (x_a + x_b)/2
            f_c, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(x_c)

            if (xd_H != maximum(abs.(f_c)))
                xd_H_state = "declining"
            else
                xd_H_state = "stalled"
            end

            xd_H = maximum(abs.(f_c))

            fc_samesign_a = abs.(sign.(f_c)+sign.(f_a))/2
            fc_samesign_b = abs.(sign.(f_c)+sign.(f_b))/2
            
            x_a = x_c.*fc_samesign_a + x_a.*(ones(params.I,1)-fc_samesign_a)
            x_b = x_c.*fc_samesign_b + x_b.*(ones(params.I,1)-fc_samesign_b)

            iter_H = iter_H + 1

            wH_i = x_c

        end

        if (xd_H <= xd_tol_H && iter_H <= iter_H_max)
            labor_market_cleared = 1
        else
            labor_market_cleared = 0
        end

        if (iter_H == iter_H_max)
            xd_H_state = "maxed-out"
        else
        end


        #### Evaluate at equilibrium wages
        xd_H_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, A_i_c, sumexp_C, pi_cC, pi_N, pi_C, pi_c, L_i_c, QS_i_c = excessdemand_H(wH_i)
        xd_H = maximum(abs.(xd_H_i))


        #### Supply elasticities
        dQdP_aux = (params.theta_lambda-ones(params.I,1)).*(ones(params.I,params.C)-pi_cC) + (params.theta-ones(params.I,1)).*pi_cC.*(ones(params.I,1)-pi_C)
        dQdP = (dQdP_aux./params.gamma_L_i) + (ones(params.I,params.C)./params.gamma_L_i)-ones(params.I,params.C)
        dPdQ= dQdP.^(-1)


        ################ Intermediaries

        #Marginal product
        mp_ij_c = reshape(params.mp,params.I*params.C)*ones(1,params.J)

        #Markdown
        if params_policy.conduct=="pc"
            mu_ij_c = ones(params.I*params.C,params.J)
        else
            s_ij_c = s_ij_c./s_ij_c
            N_ij_c = reshape(params.N_firms,params.I*params.C)*ones(1,params.J)
            mu_ij_c = ( ones(params.I*params.C,params.J) + reshape(dPdQ,params.I*params.C,1).*(s_ij_c./N_ij_c) ).^(-1)
        end


        ################ Demand

        #County-destination-commodity prices
        p_ij_c_pretax = (reshape(p_i_c_pretax,params.I*params.C,1).*repeat(params.tau_ij,params.C,1))./(mu_ij_c.*mp_ij_c)
        p_ij_c = p_ij_c_pretax + params_policy.tax_ij_c

        #Nation-destination-commodity price index
        aux_n = sum(repeat(params.a_ij_c.*(p_ij_c.^(ones(1,params.J)-params.eta_l')),1,params.N*params.C).*kronecker( Diagonal(ones(2)),kronecker(params.indicator_in,ones(1,params.J))),dims=1).^((ones(1,params.J*params.N*params.C)-kronecker(ones(1,params.C+params.N),params.eta_l')).^(-1))
        P_nj_c = reshape(aux_n,params.J,params.N*params.C)'

        #Destination-commodity price index
        aux_j1 = params.a_nj_c.*(P_nj_c.^(ones(1,params.J)-params.eta_m'))
        aux_j = hcat(aux_j1[1:params.N,:],aux_j1[params.N+1:2*params.N,:]) #horizontally concatenate C times
        P_j_c = (reshape(sum(aux_j, dims=1),params.J,params.C).^((ones(params.J,1)-params.eta_m).^(-1)))'
        
        #Destination price index
        P_j = sum(params.a_j_c.*(P_j_c.^(ones(1,params.J)-params.eta_u')),dims=1).^((ones(1,params.J)-params.eta_u').^(-1))

        #County-level demand
        aux_ij_c = params.a_ij_c.*(p_ij_c.^(-params.eta_l'))
        aux_nj_c = (params.a_nj_c.*(P_nj_c.^(params.eta_l'-params.eta_m')).*kronecker(params.a_j_c.*(P_j_c.^(params.eta_m'-params.eta_u')),ones(params.N,1))).*(P_j.^(params.eta_u'-ones(1,params.J))).*params.X_j'

        C_ij_c = aux_ij_c.*(kronecker(Diagonal(ones(params.C)),params.indicator_in)*aux_nj_c)
        Q_ij_c = C_ij_c.*repeat(params.tau_ij, params.C)
        QD_i_c = reshape(sum(Q_ij_c, dims=2),params.I,params.C)
        
        #Source shares
        s_ij_c = Q_ij_c./sum(Q_ij_c, dims=2)

        #Excess demand
        ED_ic = reshape(QD_i_c - QS_i_c,params.I*params.C)
        QDQS_den_max = maximum( hcat(reshape(QD_i_c,params.I*params.C),reshape(QS_i_c,params.I*params.C)), dims=2)
        QDQS_den_min = minimum( hcat(reshape(QD_i_c,params.I*params.C),reshape(QS_i_c,params.I*params.C)), dims=2)

        xd_ic_max = ED_ic./QDQS_den_max
        xd_ic_min = ED_ic./QDQS_den_min

        return (xd_ic_max, xd_ic_min, pi_cC, pi_C, pi_c, pi_N, p_i_c_pretax, p_i_c, p_ij_c_pretax, p_ij_c, QD_i_c, QS_i_c, C_ij_c, Q_ij_c, s_ij_c, P_j_c, sumexp_C, dQdP, dPdQ, mu_ij_c, labor_market_cleared, iter_H, wH_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, xd_H, xd_H_i, xd_H_state)

    end


else
end

####################### SOLVE MODEL FOR EQUILIBRIUM PRICES

#Solver settings
step_size = params_solver.step_size
max_iterations = params_solver.max_iterations
start_price_multiplier = params_solver.start_price_multiplier

#Starting values
if params_policy.policy== "lf"
    p_ic = start_price_multiplier*ones(params.I*params.C,1)
    s_ij_c_0 = ones(params.I*params.C,params.J)./params.J
else
    p_ic = start_p_ic
    s_ij_c_0 = start_s_ij_c_0
end
xd_ic_max, xd_ic_min, pi_cC, pi_C, pi_c, pi_N, p_i_c_pretax, p_i_c, p_ij_c_pretax, p_ij_c, QD_i_c, QS_i_c, C_ij_c, Q_ij_c, s_ij_c, P_j_c, sumexp_C, dQdP, dPdQ, mu_ij_c, labor_market_cleared, iter_H, wH_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, xd_H, xd_H_i, xd_H_state = excessdemand(p_ic,s_ij_c_0)

xd_tol = params_solver.xd_tol
xd_max = xd_tol + 1
xd_min = xd_tol + 1
n_iterations = 0
while (xd_min > xd_tol) && (n_iterations<max_iterations)

    global xd_max, xd_min, n_iterations, p_ic, s_ij_d
    global xd_ic_max, xd_ic_min, pi_cC, pi_C, pi_c, pi_N, p_i_c_pretax, p_i_c, p_ij_c_pretax, p_ij_c, QD_i_c, QS_i_c, C_ij_c, Q_ij_c, s_ij_c, P_j_c, sumexp_C, dQdP, dPdQ, mu_ij_c, labor_market_cleared, iter_H, wH_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, xd_H, xd_H_i, xd_H_state

    s_ij_c_inn = s_ij_c
    xd_ic_max, xd_ic_min, pi_cC, pi_C, pi_c, pi_N, p_i_c_pretax, p_i_c, p_ij_c_pretax, p_ij_c, QD_i_c, QS_i_c, C_ij_c, Q_ij_c, s_ij_c, P_j_c, sumexp_C, dQdP, dPdQ, mu_ij_c, labor_market_cleared, iter_H, wH_i, wH_i_c, r_i_c, rtilde_i_c, zeta_i_c, xd_H, xd_H_i, xd_H_state = excessdemand(p_ic,s_ij_c_0)    
    p_ic = p_ic + xd_ic_max.*step_size

    s_ij_d = maximum(s_ij_c - s_ij_c_inn)
    xd_max = abs(maximum(xd_ic_max))
    xd_min = abs(maximum(xd_ic_min))
    n_iterations = n_iterations + 1

end

if (n_iterations < max_iterations)
    if (labor_market_cleared==1)
        print("Case: ", params.mobility_spec, " + " , params_policy.conduct, " + ", params_policy.policy, ". Success in output and labor market. ")
    else
        print("Case: ", params.mobility_spec, " + " , params_policy.conduct, " + ", params_policy.policy, ". Success in output market, failure in labor market (", xd_H_state,"). ")
    end
else
    if (labor_market_cleared==1)
        print("Case: ", params.mobility_spec, " + " , params_policy.conduct, " + ", params_policy.policy, ". Failure in output market, success in labor market. ")
    else
        print("Case: ", params.mobility_spec, " + " , params_policy.conduct, " + ", params_policy.policy, ". Failure in both markets (labor market ", xd_H_state,"). ")
    end
    hcat(reshape(xd_ic_min, params.I, params.C), xd_H_i)
end



####################### EVALUATE MODEL AT EQUILIBRIUM PRICES

#Define output folder
if params.mobility_spec=="baseline"    
    results_folder=string("results_",params_policy.conduct,"/",params_policy.policy)
else
    results_folder=string("results_",params_policy.conduct,"/",params_policy.policy,"/mobility")
end

#Main variables
L = hcat(pi_c,pi_N).*params.L_i
Q = rtilde_i_c.*(pi_cC.^(-params.theta_lambda.^(-1))).*(pi_C.^(-params.theta.^(-1))).*L[:,1:end-1]
C = hcat(sum(C_ij_c[1:params.I,:],dims=1)', sum(C_ij_c[params.I+1:2*params.I,:],dims=1)')
mu_i_c = reshape(mu_ij_c[1:end,1], params.I, params.C)

#Save market outcomes
writedlm(string(results_folder,"/L.csv"),L,',')
writedlm(string(results_folder,"/Q.csv"),Q,',')
writedlm(string(results_folder,"/p_i_c.csv"),p_i_c,',')
writedlm(string(results_folder,"/p_i_c_pretax.csv"),p_i_c_pretax,',')
writedlm(string(results_folder,"/wH_i.csv"),wH_i,',')
writedlm(string(results_folder,"/C.csv"),C,',')
writedlm(string(results_folder,"/C_ij.csv"),C_ij_c,',')
writedlm(string(results_folder,"/Q_ij.csv"),Q_ij_c,',')
writedlm(string(results_folder,"/p_ij_c.csv"),p_ij_c,',')
writedlm(string(results_folder,"/s_ij_c.csv"),s_ij_c,',')
writedlm(string(results_folder,"/p_ij_c_pretax.csv"),p_ij_c_pretax,',')
writedlm(string(results_folder,"/mu_i_c.csv"),mu_i_c,',')
writedlm(string(results_folder,"/zeta_i_c.csv"),zeta_i_c,',')
writedlm(string(results_folder,"/dQdP.csv"),dQdP,',')

#Save policy parameters
writedlm(string(results_folder,"/subsidy_i.csv"), params_policy.subsidy_i,',')
writedlm(string(results_folder,"/tax_i.csv"), params_policy.tax_i,',')
writedlm(string(results_folder,"/tax_ij.csv"), params_policy.tax_ij_c,',')

